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Self-Consistent Quasi-Particle RPA (SCQRPA) is for the first time applied to a more level pairing 
case. Various filling situations and values for the coupling constant are considered. Very encouraging 
results in comparison with the exact solution of the model are obtained. The nature of the low lying 
mode in SCQRPA is identified. The strong reduction of the number fluctuation in SCQRPA vs BCS 
is pointed out. The transition from superfluidity to the normal fluid case is carefully investigated. 
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One of the most spectacular quantum phenomena is the transition to the supraconducting or superfiuid state in 
interacting Fermi systems. This happens e.g. in metals, liquid 3 He, neutron stars, in finite nuclei, and it is actively 
searched for in systems of magnetically trapped atomic Fermions. In most of these systems the canonical mean field 
approach of Bardeen, Cooper, and Schricffer (BCS) with a couple of adjustable parameters works astonishingly well. 
However, in recent years there have been increasing attempts to describe the pairing phenomenon on completely 
£f~) \ microscopic grounds. To our knowledge these attempts have mostly been carried out for nuclear systems. This stems 
on the one hand from the fact that phenomenological NN forces are on the market which very well describe the 
nucleon-nucleon phase shifts in all channels and in a wide range of energies. On the other hand the physics of neutron 
stars makes quantitative predictions of the pairing phenomenon in neutron matter indispensable, since superfluidity 
of neutrons stars manifests itself only quite indirectly through e.g. the phenomenon of neutron star glitches. The 
microscopic approaches to pairing, starting from a bare two body interaction, are not very numerous. The simplest 
one is based on BCS theory, using however, in the gap equation the bare force and for the single particle dispersion 
the one given by Bruckner theory. In this way one obtains e.g. gap values in the channel for neutron- neutron 
pairing which in infinite matter, as a function of the Fermi momentum Uf , have a typical bell shaped form roughly 
dropping to zero around kp = 1.3/m _1 and culminating at kp = 0.8/m _1 to values of A = 2.5 — S.OMeV for 
£h ■ neutron and nuclear matter respectively. This rather elementary approach has been extended in the past in various 
ways. The most ambitious procedure is probably the so-called correlated basis function approach [Q. However, 
more recently self-consistent T-matrix approaches and extended Bruckner theories with rearrangement terms have 
achieved a remarkable degree of sophistication j| . The screening of the interaction was treated to lowest order in the 
density, resuming the RPA bubbles, in introducing self-consistent Landau parameters ||. The outcome of all these 
investigations inevitably leads to a quite substantial reduction of pairing in neutron matter but also in symmetric 
nuclear matter. The global reduction generally attains important values and often reaches factors close to three. Such 
small values of the gap in infinite matter, however, pose a problem. Employing the Local Density Approximation 
(LDA) to estimate from the infinite matter results the gap in finite nuclei Q , one reaches with the simplified approach 
described above using the bare NN force quite reasonable gap values for finite nuclei. Interestingly in the gap equation 
quite similar results are obtained with the Gogny D1S force using the same procedure. However, with such strongly 
reduced gaps from the more sophisticated approaches mentioned above, one obtains much too small gaps in finite 
nuclei. Of course, this reasoning may be completely erroneous and the situation in finite nuclei may be very different 
from infinite matter. Nevertheless we find the above argumentation intriguing. On the other hand we know that 
pairing is an extraordinarily subtle process and employing theories which are in one or the other way uncontrolled 
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may turn out to be a hazardous enterprise. In such a situation it is probably wise to investigate the problem from 
different sides using a variety of approaches. 

In the past we have made very positive experience with an extension of RPA theory which we called Self-Consistent 
RPA (SCRPA) [||.|. For instance, in a recent work this theory has been applied to the exactly solvable many level 
pairing model in the pre-critical regime and very good agreement with the exact results for ground-state energy and 
the low lying part of the spectrum was found [@, ||] . This success has encouraged us to develop the SCRPA formalism 
also for the fully developed superfluid regime. This is a not completely trivial extension of the SCRPA and we here 
apply it for the first time to the two level pairing model. As we will see the theory also gives very promising results 
in the superfluid phase. Since the Self-Consistent Quasi-Particle RPA (SCQRPA), as in general the SCRPA theory, 
can be derived from a variational principle, which turns out to be very close to a Raleigh-Ritz variational theory, we 
believe that SCQRPA is a non perturbative approach going in a certain systematic way beyond the mean field BCS 
theory, including in a self-consistent way correlations and quantum fluctuations. It is our believe that this microscopic 
approach can ultimately be used to calculate pairing properties of realistic Fermi systems starting from the bare force. 

It should be mentioned that extensions of RPA theory, based on the Equation of Motion (EOM) method, have 
by now a quite long history. They, to a great deal, have been developed in nuclear physics. It started out with the 
work of Hara who included the ground-state correlation in the Fermion occupation numbers || . More systematic was 
the consequent work by Rowe and co-workers (see the review by D. J. Rowe fl(i|| ). The same theory was developed 
using the Green's function method by one of the present authors || . Independently the method was also proposed by 
Zimmermann and G. Ropke plus coworkers using a graphical construction pj]] . These authors named their method 
Cluster-Hartree-Fock (CHF) and it is equivalent to Self-Consistent RPA (SCRPA). The latter approach has recently 
been further developed by Dukelsky and Schuck in a series of papers [0, |, ||, g |f, |||. However, also other authors 
contributed actively to the subject A number of remarkable results have been obtained with SCRPA in non 

trivial models where comparison with exact solutions was possible |l2| . For instance for the exactly solvable many 
level pairing model of Richardson |l7j] SCRPA provides very accurate results for the ground state and the low lying 
part of the spectrum (?], || 



In detail our paper is organized as follows : in section II the two level pairing model is introduced, in section 



III 



the SCQRPA formalism is presented, in section ^ numerical results are given and detailed discussions are presented 
Comparison with other recent works is made in section [v|, in section VI the quest ion o f the second constraint on the 



particle number variance is invoked and applied to the Seniority model. In section VII , we will summarize the results 



and draw some conclusions. Finally, some useful mathematical relations and a second method for the calculation of 
occupation numbers are given in the Appendices. 



II. THE MODEL 



The two-level pairing model is an exactly solvable model extensively employed in nuclear physics to test many-body 
approximations. It was first used to test the pp-RPA and its ability to describe ground-state correlations and 
vibrations in the normal phase as well as in the superfluid phase. The model is composed of two levels with equal 
degeneracy 2ft = 2J+ 1 (J is the spin of each level) and an single-particle energy splitting e. The pairing Hamiltonian 
in this model space is 

H = \ E 3$i - & E A ) A r > J = ±! (!) 

3 33' 

where j takes the values 1 for the upper level and —1 for the lower level. Nj and Aj are the number and monopole 
pair operators of the level j, respectively, 

t 1 ° t t 

A j = 7y|| E a )rn a )rh ( 2 ) 
* rn — 1 

and 

n 

N 3 = E ( a jm a jm + Ojm a jm)- (3) 
m— 1 

where aj m creates a particle in the level j with spin projection m and ajm = (— 1) J ~ Tn <Zj_ m . The operators obey the 
following commutations relations, 

^■,41 = ^(1-1). 
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[Nj,A f ] = -SjjAAj, (4) 

thus, they define an SU(2) algebra for each level and the two level model satisfies an SU(2) x SU(2) algebra. 
For a system not at half filling, the normalized states in the Hilbert subspace of the monopole pairs are 



where f2 = f2 leads to the half filling case, i.e. the lower level is filled for 5 = 0. The matrix Hamiltonian is tridiagonal 
of dimension + 1, with matrix elements 

h ntn = (n\H\n) = e(2n - fi) - G(2n& - 2n 2 + (ifl - Q 2 + fi), (6) 



h n -i,n = (n- l\H\n) = -G<Jn(Cl - {n - l))(f2 - Cl + n)(fi -n + 1) (7) 
where, n is the number of pairs in the upper level and the number of particle is given by N — 2f2. 

III. SELF-CONSISTENT QRPA 

In a recent work || the SCRPA has been applied with very good success to the picket fence model in the non 
superfluid phase. The extension to the superfluid phase is slightly delicate and we here limit ourselves to the two 
level model, however considering arbitrary degeneracies and fillings of the levels. The objective in this section is to 
establish the equations for the Self-Consistent Quasi-Particle RPA (SCQRPA). A first application of SCQRPA has 
been performed in jl4j for the case of the seniority model (one- level pairing model). We will again later come back to 
this model. Here we want to consider the two level pairing model with arbitrary filling and coupling strength in the 
SCQRPA approach which already more or less shows the full complexity of more realistic many level problems. As a 
first step we have to transform the constrained Hamiltonian 

H' = H - fj,N, (8) 

where N is the full particle number operator, to quasi-particle operators 

°in) = ( u i - v A( a U) (9) 

OijfnJ V u i ) V".'-. / 

^A = ( u i. v A( a }A do) 



jm J \ v j u j J V a jfh 



with 



We define new quasi spin operators as 



u3+v? = l, j = ±l. (11) 

m>0 

and the quasi-particle number operator in the level j is given by, 

Nq,j = ( a jm a jm + ajmajm)- (13) 

m>0 

The quasi-particle operators obey the following commutations relations, 



[N qd ,Pj,] = 8 jr 2P], 
[N q , 3 ,Pf] = -S jr 2Pj. 

Then the Hamiltonian in the quasi-particle basis can be written as 

H' = Hqq + H' n + H 20 + H 22 + H' 31 + H' 4Q + H' n _ n 

where 



H' 00 - h , 

H' n = fciJV, i i + ft_iJV,,_i, 

H20 = ^ 2 (P 1 t +Pi) + /i- 2 (Pl 1 +P-i), 

H' 22 = ^Pi+Zi-s-PliP-i+^^P-i+PliPi), 

H' 31 = h B (P}N qtl + N q ,iPi) + hsiP^Ng,-! + N q ,-iP-i) 

+ heiP^Ng,-! + N g ,-iPi) + h-eiP^Ng,! + N q ,iP-i), 

H' i0 = h 7 {P\Pl + Pi PO + ft_ 7 (Pi 1 Pl 1 + P_iP_i) + h 8 (PlP^ + P_ 



h Q = {e-2pt)Vtvl-gn{Q.u\vl + vl)-{e + 2ii)nv\-gn{nu\v\+v i _ 1 ) 
— 2gQ 2 uiViii-iV-i, 

hi = (- - v)(u\ - vj) + gQ,(2u\v\ + -±) + 2gtt 2 uiViii-iV-i 7 

e v 4 
h-i = -(-+ii)(v 2 _ 1 -v\)+gn(2u 2 _ 1 v 2 _ 1 + -^-) + 2gn 2 U-iv^iUiVi, 

h 2 — VT}uiVi(e — 2fi) — gQ^uiVi(u 2 — v\)VVl + ^y^ 1 | ~ g^VTiu-iV-i( 

m-i«-i(m 2 _ 1 — v^^VCt H J 



— gQ,VciuiVi(u 2 _ 1 — v 2 _ j), 

/i 3 = — + Ui), 

/i_ 3 = -gfi^! + vti), 

hi = —gQ(ulu 2 L 1 +v 2 v 2 _ 1 ), 

h 5 = gVUuivi(ul - v\), 

h-5 = gVUu-iV-i(u 2 _ 1 — vt^), 

h 6 = gVUu-iV-i(ul - v\), 

h-6 = gVUuiVi(u 2 _ 1 — w^i), 

h 7 = gSlujy 2 , 

h-7 = gClu_ 1 v_ 1 , 

h$ = gfl(u 2 1 v 2 _ 1 +u 2 _ 1 v 2 ), 

hg — —g£lu\v 2 , 

h-9 = —gCtu_ 1 v_ 1 , 

hio = —2gSluiViU-iV-i. 

Also, in this basis the full particle number operator is given by, 

n = J2^, i = ±i 

3 

where, 

Nj = (u 2 - v?)N qJ + 2Q,v 2 + 2ujVjVQ,{P] + Pj). 



5 



(27) 



The RPA excited states are, as usual, obtained as 

\u)=Qt\RPA), 

where \RPA) is the correlated RPA ground-state defined via the vacuum condition 

Q V \RPA) = 0. (28) 

In terms of the generators of the Hamiltonian Ng,j, Pj and Pj, for the most general QRPA excitation operator, which 
can be viewed as a Bogoliubov transformation of Fermion pair operators 1 , we can write down the following expression 



Ql = Xj,»P} - YfrPj, " = 1, 2, 
i=±i 



where we introduced the following notation, 



1 



=, i = ±i, 



(29) 



(30) 



guaranteeing that the RPA excited state ( |27j ) is normalized, i.e. {v\v') = 8 v y . The RPA amplitudes Xj M and Yj tl/ 
in (29) shall obey the following orthogonality relations, 



E ^-^ = i." = i.2 

j=±i 

X-i^X-i .2 + Xi t iXi t 2 — Y-\_\Y-\_2 — 5^i,iii,2 = 

X\ 2*1,1 + -X-l,2*-l,l — -Xi, 1*1.2 — -X_l,l*-1,2 = o 



and the closure relations, 



E x 

v=l,2 



Y 2 

].y 3,v 



1, j = ±l, 



(31) 



A-i^Aij + X—xpXift — Y—\ t \Yxj. — *-i,2*i,2 — 0, 



with which one can invert relation (E9 











Pi 


■( 









1,2 — A_ 


-1,1*1,1 


— A-1,2 


Ai,i 


Xl,2 


Yi.i 


n,2 


*-M 


-X- 1,2 


Y-i,i 


^-1,2 


*M 


^1,2 


x i,i 


Ai,2 


y-i.i 


^-1,2 


A_i.i 


-X^-1,2 



Ql 
\QlJ 



(32) 



(33) 



In analogy to Baranger ]l9| we obtain the SCQRPA equations in minimizing the following mean excitation energy 

{[Qu,[H',Qi]}) 



[Qu, Ql]) 



(34) 



with respect to the RPA amplitudes Xj }V and Yj >v . The minimization leads straightforwardly to the following eigen- 
value problem 



1 We can not include the Hermitian pieces N q ,j in (^), since this leads to non-normalizable eigenstates as in the case of Goldstone modes. 
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/ 




Ai, 2 




Bl,2 \ 


( x h» 


A 2 ,i 


A 2 ,2 










-B 1A 


—Pi, 2 




:¥;) 




V 


—-82,1 


— B2,2 


-A 2A 







( *l,u 



(35) 



where, 



A1.1 = {[P,,[H',PI]}), A 1 , 2 = ([P 1 ,[H',P^}]) 

A 2A = ([P. u [H',Pt]}), A 2t2 = ([P_ 1 ,[H',pi_ 1 }]) 

Bi,i = -([Pi.t^A]]), B 1 ,2 = -([P l ,[H',P-i]]) 

B 2A = -([P-xJ^Pi]]), B 2 , 2 = -([P. 1 ,[H\P- 1 ]}), (36) 

and the (...) stands for the expectation values in the RPA vacuum defined by (p8|). Explicitly, the RPA matrix 
elements are given by 



r 9 / P t P \ 1 _ gjjW) 4. 1 9 1 r 



1 1 Mm) 1 Mm) / iiM 

-, M,i) + M.-i) , (N q ,iN q .-i) / P pt \ 



A h2 = h 4 " : +4/t iu 

1 M,i) + M,-i) M,iAV-i) /p p f> 

A a>1 = ^ g + ,^ + 4ft 10 - <P - lPl> 



9 /pt p \ I 2(jV,,_ 1 ) Mg,-i/ n „ , 



n 



I 1 _ Mg,-l) 1 _ ("g.-l) J 

x a 1 n 



f2 

f 1 

2/i 7 



B M = -o 77T^ ^(PiPi) + /i4(P-iPi) + hsiP^Pi) } + 8h 



1 (ptpo + (p 1 pI ) , 1 - ^ + ^ 



1 _ (N q .i) 1 _ M.i) 

1 q 1 n 



n 



-1 M<m) + Ms,-i) , Mg.i^g.-i) /p p x 

P 2 ,! = ft, *" n . + ~ 1 ^ +4h 10 - {PlP - 1 ' 



tt- = -J i Jh- 3 {P-iP-i) + h 4 (P 1 P- 1 ) + h s (PlP- 1 )\ + 8h- 9 {J ' '' r ' ' 

1 _ (iVg.-l) I 

- 1 - r~i 



Qi (^,-i> I ' ' 1 'J i W.-i) 

1 n 



n -1 _ (iV g ,-i) 1 _ <AT g ,-i) 

1 n 



Using d33| ) and the condition ( p^ ) the expectation values of type (P-Pj), {PjPj), (P-P-) and (PjPj) are readily 
expressed by the RPA amplitudes and „ (these calculations are detailed in Appendix [a]) 
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(PlPi 

(PiP} 
(P-iPU 
(PM 
(Pi Pi 



(P'-iPU 



(P-1P-1 
(PiP-i 



(Ptp_ 



(PxPU 



(n,!+n, 2 )(i-%^) ; 



= (X\ 1+ X\ 2 ){1- 



(x 1 , 1 y 1 , 1 + x 1 , 2 r 1 , 2 )(i- 



(x_ 1>1 y_ 1 , 1 + x_ l!2 y_ 1>2 ) (i - ^g^), 



= (A-_i i iy_ 1 ,i + A-_ 1 , 2 y_i l2 )(i- 



8,-1, 



(X 1A Y_ M + X!, 2 Y_ 1>2 W (l - 



= (*_ M Y M + A^Y^W (l - 



(y 1|1 y_ lil + y 1 , a y_ 1 , a )J(i-<M)( 1 



(x M A_ ia + x 1>2 x_ 1)2 ) J(i- { -^A) (i - . 



(38) 



Before we discuss how to express the expectation values (N q j), (JV? j) and (N q jN q ji) as functions of the amplitudes 
Xj v and Yj v , we want to give the equations for the determination of the Bogoliubov amplitudes Uj, Vj of equations (|^) 
and ( |Io| ) . As usual they are determined from the minimization of the ground-state energy |L3| p0[ 

d(H') d(H') dvj 



du-i 



dvj duj 



([H',P]]) = 0, j = ±l 
([H',Qt}) = 0, i/ = l,2 



(39) 



with 



(H 1 ) = ho + + ft_i(^ 9l _i) + ^(P/Pi) + ^(PliP-i) + ^((P/P-i) + (PlxPi)) 

+ ^((Pi 1 ^) + (PiPi)) + h- 7 ((P^-i) + (P-iP-x)) + M(PlPU) + (P-iPi)) 

+ h^NlJ + /i_ 9 <A\,_i) + h w {N qil N qi _j). (40) 

It should be mentioned that when ([59]) is evaluated with a BCS ground-state then (|3^) leads to the usual BCS 
equations. However, here we use the correlated RPA ground-state and then the mean field equations couple back to 
the RPA amplitudes Xj v and Y^ v . Explicitly these equations lead to 



2^u jVj + A,(«2 - u)) = 0, j = ±1 
which together with ([ll]) can be written as 



J - 



(41) 



(42) 



with the standard solution 
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^ = V 1 - -m L= ^ ) (43) 



from where follows the gap equation, 



Ai = !'-j": r j = \ 9ij i 3 =, h 3 = ±1 (44) 



where the renormalizcd single-particle energies are 



= (if - 9$ + — |— (« 2 -, - ^)«^> + (PjP-i)) - M, i = ±1 (45) 



1 - 

and the renormalizcd interaction is given by 



n 



with, 



9i,i 



g = ( f; 1 f- 1 ) , (46) 
5-1,1 9-1,-1 



g n- ^ i) {2((p 1 p 1 ) + (ptpx)) + (iv ?!l ) - 

5i,-i = .9^ 



1- n 



n 



i 



Q 

(N ■,) - iMdl 
9*1 -g- 



1 



f2 

2 



ff n - ^ {2((p! 1 p! 1 ) + (PliP-i)) + - ^ (47) 



l 



We see that the mean field equations have exactly the same mathematical structure as in the BCS case, however, 
with renormalized vertices and single-particle energies involving the RPA amplitudes. We, therefore, explicitly see 
that the mean field equations are coupled to the quantum fluctuations. 

Let us now come to the elaboration of the quasi-particlcs occupation numbers and their variances. The determination 
of those quantities is one of the difficulties in the SCQRPA approach ^ [n| ^o). However, this problem has found an 
elegant solution in the early works of j2lj (see also ||22) ) . In the same way, we derived expressions of the quasi-particles 
occupation numbers and their variances as expansions in the operators P- and Pj up to any order in a systematic 
way. The detailed derivation is given in the Appendix |^. We here present a different method which shows some 
interesting aspects and will lead to the same result. Using the bosonic representation of the quasi-spin operators of 
our model, we can write 





= 2P]P,-, 


P} 


= B)(l- 


Pj 


= (P}Y = 



Q 1 3 
1 



Bj (48) 



where, one can show that these operators in this representation always obey to the commutation rules of angular 
momentum (pT|). We also can invert this relation, and we obtain 



B\ = PU 1 - ^B\B 



1 



3 3 \ <> ■< '• 

Bj = (l-ip)s^ 2 Pj . (49) 
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With (E9J) N q .j can be expressed as 



K,3 = ZflJ^i, 

-1 



= 2pni--B^ />. 



- ^-i^-) lp - (50) 

Therefore, we obtained a recursive relation for N q ,j, and with it we can derive an expansion for N q ,-. By successive 
replacement of iV g j in the rhs of j5C|), one finds the following expansion, 

+ 2 +2 ~ f^p]Pj - Ngj + n\ n „ 



' ttt' • •••• ( 51 ) 



It should be noted that the first term in ( |51| ) becomes already exact for J = 1/2 and including the second term it is 
also exact for J = 3/2, etc. 



For JV| j , we can use the Casimir relation, 



np i P . + n + i^ = o. (52) 



It is equivalent to use the expansion of N q • obtained as the square of iV gi . 



^=4ptp J + ^±^^ + .... (53) 

In the same way, we use (|5l| ) to obtain an expansion for N q ^N q ^i, but it is sufficient to use the term of the first 
order of this expansion, to obtain 

N^iNg,-! = tPlPiP^P-i + ... (54) 

In principle the expansion ( |5l| ) can be pushed to higher order, however, it quickly becomes quite cumbersome and in 
practice we always will stop at second order. In any case the expansion is finite with maximal J + 1/2 terms. It is 
natural that such an expansion exists since there is a duality between the pair of operators f?j , Bj <^4> P] , Pj . There 
is the choice either to bosonize the problem then everything is expressed in terms of £>J and Bj operators. Or one 
stays with the Fermion pair operators and everything is expressed in terms of P^ and Pj. In [p3| the former route 
was chosen, here we choose the latter one. One should mention that a truncation of the series (pW) also entails some 
violation of the Pauli principle but one may notice that the series is very fast converging and that already the lowest 
order correctly contains two limits : J = 1/2 as already mentioned and J — > oo, since then Pj — > Bj and the lowest 
order is also correct see (fl8|). With these remarks in mind we go ahead. By the inversion of the QRPA excitation 
operator Qj,, the expectation values of these expressions are immediately given in terms of the RPA amplitudes Xj u 
and 5^,1/, as one can see in appendix ^ where we give some details concerning the calculation of expectations values 
of these expressions in the RPA ground-state. 

Our system of SCQRPA equations is now fully closed and we can proceed to its solution. Before let us, however, 
shortly come back to the limit of standard QRPA. This we will do for the symmetric case i.e. N — 20. This case 
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is obtained in evaluating all expectation values in all interaction kernels with the BCS ground-state or else putting 
Yj v =0 and ^ y X? = 1 for j = ±1. The matrix elements are then 

A 2 A 2 





= ^2,2 


Al, 2 


= A 2 ,l 


Bi,i 


= B2,2 


-Bl,2 


= -02,1 



A 2 
2gtt' 

A 2 A 2 
2gQ 2 + IgU 1 
A 2 



where, the gap equation in the BCS theory leads to the solution in the symmetric case 



2 



together with, 



A = y ff 2 r! 2 -^, (56) 



1 



it; = vt t = -I i 



2 2 
V x = U_t = 



2 \~ ' 2 g nJ' 



M = -! (57) 



2 V 2 5^/ 

2 

where £ is defined as £ = 2ef2/ (2f2 — 1). For the positive eigenvalues of the RPA matrix, we obtain 

Sl? RPA = 0, (58a) 
nr PA = ^A 2 -^. (58b) 

As usual the other two eigenvalues are 

_qQRPA with v = 12 . These results are well known [fl8| , |24|| . We have 
repeated them here for completeness and stressing the point that in QRPA, because of the spontaneously broken 
particle number symmetry, one obtains a Goldstone mode flf RP = 0. We again would like, to stress the point that 
this is the case only if we evaluate (55) with the solution Uj, Vj given by the mean field equations ( ffl| ) which for 



Y^ V X? — 1, Yj> = reduce to the usual BCS equations. We explicitly showed it here for the symmetric case but 
the same scenario holds true for cases away from half filling. 



IV. RESULTS AND DISCUSSIONS 



We first recall that the phase transition point in BCS theory for the two level pairing model is produced at 
g c = e/(2fi — 1), where e is the single-particle energy splitting and f2 is the pair degeneracy of each level. In the 
following, the graphs are plotted, as usual, as function of the variable V — gft/2e, and refer to the case with level spin 
J = 11/2, i.e. = 6 and single-particle energy e = 2 (in arbitrary units). This latter value for J has been chosen for 
easier comparison with the results of ||| which will be given in section 

Let us first discuss the case with N = 12, i.e. the lower level is filled in the absence of correlations. We call 
this the half-filled or symmetric case. In Fig.^ we show in the upper panel the excitation energies. Let us consider 
the well known scenario of the standard RPA. Before the phase transition to the superfluid phase we work with 
the unconstrained Hamiltonian. One obtains two eigenvalues with the interpretation of differences of ground-state 
energies, differing by two units in mass 2/i ± = ±(Eq ±2 — E^). They are evidently related to the chemical potential 
and in standard particle-particle RPA (pp-RPA) they are given by 

= 2/i+ = -g + */g~T~e^/e + g{l - 2ft), (59a) 
fir = -2yT =g+ sfg + l + g(l - 2Q), (59b) 

where Q a and fl r correspond to the addition and removal phonons of the pp-RPA, respectively. In Fig|l] the case 
— 2/i~ is shown and we will discuss the case 2/z + separately below in Fig.pl. We see on the graph the usual result, 
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namely that —2/i~ drops to zero at the phase transition point (strictly speaking only in the large limit). After the 
phase transition point we work with the constrained Hamiltonian (|l5|) in the BCS quasi-particle representation. The 
QRPA eigenvalue (58b) is also shown in Fig.[j]. The Goldstone mode ( |58a| ) at zero energy corresponds to a rotation 
in gauge space whereas the second eigenvalue corresponds to the "/3— vibration" of the nucleus with N particles [^B] , 
This difference in interpretation is also well born out in the SCQRPA in comparison with the exact solution. We see 
that in the transition region SCRPA shows a tremendous improvement over RPA and that SCRPA follows the exact 
value of — 2/i~ even far beyond the phase transition point (as defined by BCS theory) where no RPA solution exists. 
It is also to be noticed that the sharp phase transition seen in RPA-QRPA is an artifact of the theory and that in 
reality the phase transition is completely washed out due to the finiteness of the system. The fact that the " spherical" 
SCRPA solution co-exists with the "deformed" SCQRPA solution over a wide parameter range representing different 
energy states of the system is a quite unique situation. In all other model cases where we have investigated the 
" spherical-deformed" transition the " spherical" solution ceased to converge numerically |2(J beyond a certain critical 
coupling. This, however, is no proof that the "spherical" solution does not also exist far in the deformed region 
representing physical states. It may be that in those works simply the method for the numerical solution was not 
sophisticated enough. This is a point to be investigated in the future. In the superfluid (deformed) region SCQRPA 
still is superior to QRPA but the improvement is less spectacular. This stems from the fact that the transformation to 
BCS quasi-particles effectively accounts already for some supplementary correlations in QRPA and thus the differences 
with exact and SCQRPA solutions become less important than in the non superfluid regime. A feature which is to be 
remarked in Figjl] is the fact that SCRPA and SCQRPA do not smoothly match in the transition region whereas RPA 
and QRPA have a certain continuity at the transition point. However, we see that SCRPA and SCQRPA describe 
two physically very distinct states which do not have any contact in the exact case neither and therefore it is not 
astonishing that SCRPA and SCQRPA do not join. This mismatch has as a consequence that there also exists a 
rupture in the ground-state energy as a function of interaction as is seen in the lower panel of Fig.[j]. Again SCQRPA 
results improve strongly over BCS ground-state energies in the deformed region. 

So far we have omitted the discussion of two items of the case considered in Fig.[l] which are slightly subtle. The first 
is the fact that the QRPA shows two eigenvalues : the "/3— vibration" and the Goldstone mode at zero energy ("the 
pair rotation mode"), whereas we have not shown the corresponding low energy mode of SCQRPA. We will below 
devote an extra paragraph to this issue. The second point is that we have not shown in Fig.|]the QRPA values for the 
ground-state energies. We show this separately in an enlarged scale around the transition point in Fig.|| We there see 
that QRPA overbinds in the transition region but that further to the right of the transition region QRPA values are 
closer to the exact solution than the ones from SCQRPA. This is a paradoxical result which systematically repeats 
itself for all other configurations we will consider below. However, the seemingly "better agreement" is an artifact of 
the QRPA which has already been encountered in others cases p6|. W e want to argue as follows : SCQRPA is in 
itself a well defined theory, resulting from the variational principle ( p4| ) for two body correlation functions. One also 
can consider it as a HFB approach for Fermion pairs. The Pauli principle is respected in an optimal way, since at no 
point a bosonization of Fermion pair operators is introduced and the Pauli principle is only violated in the truncation 
of (^) which is a very fast converging series. However, any approximation to the full SCQRPA scheme necessarily 
diminishes the respect of the Pauli principle what simulates more correlations than there should be. Since for the 
present model case the SCQRPA ground-state energy is systematically above the exact one (under binding), it may 
happen that, when the Pauli principle constraint is released in going from SCQRPA to QRPA, the corresponding 
gain in energy is such that, accidentally, the QRPA ground- state energy practically coincides with the exact values 
over a wide rang of parameters. We think that this is what happens in this model not only for the configuration in 
Fig.|l| but systematically for all types of degeneracies and all fillings. We will not discuss this issue for the other cases 
any more in this work. We again should mention that we have found such fortuitous coincidences already in other 
works p6| . However, in more realistic cases ones usually finds that the standard RPA strongly overbinds with respect 
to the exact values p?| . 

Let us now discuss situations where either the lower or upper levels are only partially filled. Like in the one level 
pairing case these configurations always show a non trivial BCS solution, i.e. they are always in the superfluid regime 
independent of V. Let us look at Fig.|^ with J = 11/2 and N — 8 that is the lower level partially filled for V = 0. 
In the upper panel the high lying eigenvalue of the SCQRPA equations is shown against the exact value. We see 
that there is some improvement of SCQRPA with respect to QRPA but it is not spectacular. It is similar to the 
case of Fig.[l] where in the superfluid region the improvement, for reasons already explained above, is modest. For 
the ground-state energy there is quite strong improvement over BCS theory. The QRPA result is not shown but the 
situation is the same as already explained above. The cases J = 11/2, N — 4 and N = 14 shown in Fig^-|| are 
qualitatively similar. 

Let us now come to the low lying eigenvalue of SCQRPA which in QRPA corresponds to the zero energy eigenvalue 
(Goldstone or spurious mode). In Fig.^ we show the low lying eigenvalue for the case J = 11/2 and N = 10. We 
see that this eigenvalue follows very precisely the difference 2(/i+ — fj,~) = Eq +2 + Eq~ 2 — 2Eq of the two chemical 
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potentials 2/i + = Eq +2 — E$ and 2yT = Eq — Eq~ 2 as obtained from the exact calculation. This identification 
makes indeed sense : since we are in the symmetry broken phase the SCQRPA system can not distinguish between 
N ± 2 states. For large N both 2fi + and 2/i~ tend individually to Goldstone modes but for finite N it definitely is 
reasonable to define the difference between 2/i + and 2/i~ as the low lying excitation and it is this combination which 
shows up as low lying mode in the SCQRPA calculation. This is confirmed in looking at other configurations : in 
Fig.0 we show the case J = 11/2, N = 4 and in fact we find analogous scenarios for all configurations we investigated, 
besides one : this is the symmetric case with J = 11/2, N = 12. In the Fig.|| we see that the picture is slightly 
different from the rest of the cases. This stems from the fact that in the symmetric case we have a transition from 
the superfluid to the non superfluid regime which is absent in the other partially filled cases. We also see that the 
values for 2^ + and 2/i~ are very asymmetric, 2/i + apparently taking the role of the Goldstone mode alone. Also the 
agreement of the low lying SCQRPA solution fix is slightly less good than in all other cases. 

Let us also add some remarks why in SCQRPA there is, contrary to QRPA, no exact Goldstone mode at zero energy. 
This is relatively easy to understand : in quasi-particle representation the number operator is given by (]25|), (p6[). One 
can check that in QRPA the terms a^a, if they were included, completely decouple of the QRPA equations. Therefore, 
in QRPA it is as if one had used the full particle number operator and therefore a particular solution of the QRPA 
equations is = N and with [H, N] — we get the zero eigenvalue in the Equation of Motion (EOM) approach. 
This argumentation is no longer true in SCQRPA where the terms N q .j (|2^) of the number operator contribute in 
principle to SCQRPA. However, we can not include them in the RPA operator because these are hermitian pieces 
leading to non-normalizable eigenstates. Therefore Q> = N as a particular solution only holds in QRPA but not in 
other cases such as SCQRPA. However as a benefit, we see in the preceding figures that we can identify the finite 
value of f2i with a particular rotational frequency in gauge space of the exact solution of the problem. On the other 
hand in realistic situation one can include in the RPA operator terms of the form ct\ak> for k ^ k' 
hermitian operators a\ctk have to be excluded for the reason already mentioned. These components correspond in 
an infinite system to momentum transfer zero and they are thus of zero measure. Therefore in an infinite system we 
have again full restoration of symmetry. 

Other quantities which are interesting to be calculated within the SCQRPA formalism are the chemical potentials 
directly from differences of ground-state energies. For example in Fig.p[]lO| we show /i* = ±^(Eq ±2 — Eq) where 
the individual ground-state energies are obtained directly from separate SCQRPA calculations. We see for J = 11/2 
and N = 4 and 8 that the agreement between SCQRPA results and exact values is excellent and in any case a 
strong improvement over BCS theory can be noticed. The same is true for the chemical potential fi as obtained from 
fx = 5 + fjT ) in the exact calculation. This latter which is an average chemical potential should be identified with 
the Lagrange multiplier \x used for restoring the symmetry of the good particle number (||) in BCS and SCQRPA. 
This identification is shown in each upper panel in FigJ^-|ll]. In Fig.|ll| we show the results for /i and /i* for the 
symmetric case J = 11/2 and N = 12. We see that again the same remarks as for the asymmetric cases hold true. 
However, we notice the particular situation that for fj, the exact, BCS, and SCQRPA solutions coincide exactly. This 
has to do with the specific symmetries in the half- filled case. 

It is also interesting to show the chemical potentials 2/i + and 2jjT in a symmetric way as done in pq| . This also 
gives us the occasion to study the accuracy of our approximation ( pl| ) and (^4|) for the occupation numbers. Lets 
us first of all say that we have here a quite unusual situation for SCRPA : the solution in the spherical, i.e. non 
superfluid basis, exists far into the superfluid regime. Usually in other models the solution of SCRPA in the spherical 
basis can be found up to interaction values slightly beyond the mean field transition point but here very reasonable 
values for the chemical potentials 2/i ± are obtained for all values of V as seen in Figjl^. This was also found in the 
work by Passos and al. [p8[ . It should be mentioned, however, that maintaining the spherical basis gives much less 
good results for the ground-state energy as seen in Fig.|l|. Indeed after the transition point the ground-state energy 
values deviate quite strongly from the exact results. In Fig.|l2| we calculate the expectation values (JV,-) and (NiNj) 
with the exact RPA ground-state JT|| 

\RPA) = j2(^) l (4) l (Ai 1 f- l n, (60) 

1=0 

where X, Y are the RPA amplitudes, defined with the addition (P) and removal (R) phonons of the particle-particle 
RPA and satisfying the normalization condition X 2 — Y 2 = 1. This gives the broken lines. If we calculate the same 
values from our limited expansion (|5l|) then the dotted lines are obtained. We see that beyond the transition point 
the solution becomes extremely sensitive to approximations. Indeed our approximated values deviate quite a bit from 
the ones calculated with the full wave function \RPA). 

A quantity which is particularly sensitive to the correct treatment of correlations are the occupation numbers. For 
example, for the particle number in the upper level we obtain 

(N 1 ) = (ul-vl){N q , 1 ) + 2nv 2 1 (61) 
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and the result is shown in Fig.|l3| for the superfiuid and non superfluid regimes. Once again we see that the change 
around the phase transition is not continuous. Still with SCQRPA one notices a tremendous improvement over 
standard QRPA for which the amplitudes diverge at the critical point. Indeed it is just in such quantities as occupation 
numbers where the full superiority of SCQRPA over its linearized version of QRPA is fully born out. Before finishing 
this section, we will explain how we proceeded to make the QRPA and RPA calculation of (N q ,i) in both regions 
normal and superfluid. We use the first order of the bosonic expansion of the Nq,i, i.e. the first order of the expansion 
shown in (|5l|), where it is sufficient to put P^ — P>\. Thus, with the commutation rules (14), we find 



2(Y'+Y i 



^- wmrfsiY (62) 

In linearizing this expression, we obtain 

(N qA ) = 2(Y 2 1+ Y 2 2 ). (63) 

It is interesting to detail this calculation, since it is useful to see analytically the QRPA and RPA results for the 
particle number in the upper level close the transition point. It is well known that the two excitation modes in 
the RPA method converge to zero at the transition point, then the corresponding RPA amplitudes tend to infinity, 
what explains the divergence of (Ni). In the superfluid zone, we mention that we neglected the RPA amplitudes 
corresponding to the Goldstone (spurious) mode when we make the calculation of (Nx). 

A constant concern for superfluidity or superconductivity in finite systems is that the quasi-particle transforma- 
tion ([)]) does not preserve good particle number. Even though one fixes particle number in the mean with the help of a 
Lagrange multiplier, the contamination of expectation values with components which have wrong particle number can 
be quite important. This is for instance the case for atomic nuclei. That is why, very early, one has thought of how to 
improve BCS theory with respect to particle number conservation. One quite popular approach is to project the BCS 
wave function on good particle number. An approximation to this relatively heavy scheme is the approximate particle 
number projection by Lipkin-Nogami |j29f . It is therefore interesting to investigate how much SCQRPA improves on 
the spread in particle number. We therefore will calculate 

(AA) 2 = (N 2 ) - (N) 2 (64) 



with N = Ni + N^i the particle number operator and Nj given by (|26j) within SCQRPA. The terms involving bilinear 
forms in Pj , Pj are as usual directly expressed by the RPA amplitudes and for the quasi-particle occupation number 
operators we use ([5lj)-([54"|). Then A A can be calculated and the results for various configurations are shown in 
Fig.|lHl6[ We see that the spread in particle number is strongly reduced over BCS values reaching typical factors two 



to three. We, however, see that AN even in SCQRPA acquires non vanishing sizable values. This is an expression 



that particle number is not completely restored. We will see in section VI how one eventually can improve on this. 
We also tried to evaluate AN in standard QRPA in applying a lowest order bosonization of the expression. However, 
due to the non-normalizable Goldstone mode we ran into troubles with this procedure and could not reach a definite 
conclusion on this point. 

Another interesting aspect which can be studied with our model is the question whether the pairing correlations, 
with respect to BCS, have been enhanced or weakened due to the SCQRPA correlations. To this end we define the 
following quantal expression for the correlation function 

c =^e (<^>-i^K^>y (65) 



This expression reduces to the following expression when evaluated with the BCS ground-state 

Cbcs = ^2 u]vj. (66) 

3=±1 



Multiplying (^36|) with g 2 fl 2 yields the standard BCS gap squared. We, however, refrain from multiplying (^5|) or (|6? 
with Q 2 f2 2 , since the renormalized gap from ( p4| ) is level dependent. Often also ( |65| ) is given in a non diagonal 
form [|30f| but having difficulties to express non diagonal densities with SCQRPA amplitudes we will not consider the 
non diagonal form here. We therefore evaluate ( |65| ) in three approximations : we can express (|6^) in terms of Pj , 

Pj and N q j operators and then take the expectation value with the SCQRPA ground-state. The equations (|3^), 
(pT|) and (M) then allow to express C in terms of the SCQRPA amplitudes Xj tU , Yj M . We will call this Cscqrpa- 
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TABLE I: 


Results for the ground-state 


energy (in arbitrary units) 


vs the variable V 


= gQ/2e described 


in the text. The spin 


of the levels is J = 19/2 and the number of particles is N = 20. 








V 


Exact 


SCQRPA 


BF-RPA 


QRPA 


BCS 


-0.50 


-18.55446 


-18.55410 


-18.55360 


-18.56890 


-18.00000 


-0.45 


-18.66849 


-18.66821 


-18.66784 


-18.67924 


-18.20000 


-0.40 


-18.78706 


-18.78686 


-18.78660 


-18.79474 


-18.40000 


-0.35 


-18.91072 


-18.91058 


-18.91040 


-18.91592 


-18.60000 


-0.30 


-19.04010 


-19.04001 


-19.03990 


-19.04338 


-18.80000 


-0.25 


-19.17600 


-19.17594 


-19.17588 


-19.17787 


-19.00000 


-0.20 


-19.31939 


-19.31936 


-19.31933 


-19.32031 


-19.20000 


-0.15 


-19.47153 


-19.47151 


-19.47150 


-19.47188 


-19.40000 


-0.10 


-19.63406 


-19.63405 


-19.63405 


-19.63415 


-19.60000 


-0.05 


-19.80919 


-19.80919 


-19.80919 


-19.80919 


-19.80000 


0.00 


-20.00000 


-20.00000 


-20.00000 


-20.00000 


-20.00000 


0.05 


-20.21101 


-20.21101 


-20.21101 


-20.21102 


-20.20000 


0.10 


-20.44921 


-20.44918 


-20.44917 


-20.44953 


-20.40000 


0.15 


-20.72625 


-20.72599 


-20.72593 


-20.72899 


-20.60000 


0.20 


-21.06339 


-21.06130 


-21.06100 


-21.08080 


-20.80000 


0.25 


-21.50260 


-21.48733 


-21.48640 


-21.64174 


-21.00000 


0.30 


-22.12491 


-22.03638 


-22.03620 


-22.13484 


-21.37193 


0.35 


-23.03321 


-22.80453 


-22.77899 


-22.99299 


-22.21880 


0.40 


-24.24609 


-23.99588 


-23.97001 


-24.21285 


-23.37895 


0.45 


-25.68929 


-25.42829 


-25.39821 


-25.65779 


-24.74795 


0.50 


-27.29077 


-27.01885 


-26.98390 


-27.25633 


-26.26316 



We also evaluate (|66|) in the standard BCS approximation which is ( |65|) . However, we also calculate (|66|) with Uj, 
Vj amplitudes from the renormalized BCS (r-BCS) theory, i.e. from ( |43| ) with A* solution of (fi~4|). The results are 
shown in Fig.[l7] and Fig.[lJ for N — 12 and N = 14 respectively (the case N = 10 gives exactly same results as 
N = 14). We see that r-BCS gives with respects to BCS less correlations. Eventually this suppression of pairing can 
be put into analogy with gap suppression in infinite neutron matter from renormalized theories || (see discussion 
in the introduction). However, the suppression of pairing correlation in r-BCS is misleading in our model, since, 
on the contrary, the full SCQRPA gives mostly an enhancement of pair correlations with respect to BCS. It is not 
obvious whether this conclusion can be taken over to the infinite matter case. It may, however, be indicated that the 
renormalized gap equations from screening (RPA) type correlations should be carefully treated consistently with the 
evaluation of two body correlation function before definite conclusions can be reached. 



V. COMPARISON WITH OTHER RECENT WORKS 



The two level pairing model has recently served as a testing ground for various generalizations of BCS theory. In 
spirit the work which comes closest to the present is the one of Sambataro and Dinh Dang ^3|. Instead of treating 
quasi-particle pair operators directly as we do here, they bosonize them (with a method developed in |p3[ ) and 
expand the Hamiltonian (|l|) in terms of these Bosons up to fourth order. A Bogoliubov transformation of the Boson 
operators quite analogous to our Bogoliubov transformation of Fermion pair operators (^|) is then performed and the 
corresponding non linear Hartree-Fock-Bogoliubov equation are written down. Again they are quite analogous to our 
SCQRPA equations. The coefficients of quasi-particle transformation are obtained, as usual, in minimizing the ground- 
state energy (see also our procedure (^9|)) with respect to the transformation coefficients. As in our case equations 
are obtained which couple back to the bosonic HFB, i.e. the RPA amplitudes. The coupled system of equations for 
fermionic and bosonic transformation amplitudes is then solved self-consistently. For better comparison we actually, 
on purpose, have chosen most of the configurations in J and N the same as in Q. Since in Q J = 11/2, what is 
a rather high degeneracy of the levels, the Fermion pair operators are quite collective and a bosonization certainly is 
a valid procedure. Not unexpectedly, therefore, the results of the present work are very close to the ones presented 
in p3|| . A detailed comparison shows that our results are systematically closer to the exact ones by a very small 
amount. This may be due to the fact that we never bosonize and treat the Fermion pair commutation rules exactly 
but the difference is too small for drawing any definite conclusion. In [p3| , Sambataro and al. show an explicit 
comparison of results referring to the symmetric case with J = 19/2. We also can make such a comparison for the 
ground-state energy referring to the same configuration. It is given in TABLE.[j| where we show the results for four 
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different many-body approaches : our approach (SCQRPA), approach of Sambataro and al. (BF-RPA) ||]|, standard 
Quasi-Particle RPA (QRPA) |24| and the BCS method f2q| . Of course, we recall that in our approach, we use the 
Self-Consistent particle-particle RPA in the normal fluid zone, while, in the superfluid region we use the generalized 
version of SCRPA that is SCQRPA. In order to accentuate the differences one would have to go to configurations with 
much lower degeneracies where the constraints from the Pauli principle become much more severe. For example the 
SCRPA approach has been applied to the case J = 1/2 with N = 2 in || and the exact result was recovered. It would 
be interesting to see how the approach with the bosonization pq ] performs in that case. In spite of being very similar 
in general spirit to the work in |p3[ we have solved and considered several additional problems which remained open 
in pjj| |. In first place this concerns the low lying eigenvalue of SCQRPA. No interpretation of this important root was 
given in ^3|. We, however, suppose that the results in |2^] for this state (no numerical values have been given) can 
be equally interpreted as the difference 2(/i + — as in our case. Another quantity which was not considered in |25| ] 
is the number fluctuation. Again we believe that corresponding values would be close to the ones found here. Also 
the transition to the non superfluid regime has not been treated in p3| . However, probably all these aspects will be 
quite similar in both approaches as long as a bosonization of the Fermion pair operators is valid. We think, however, 
that it does not cost much to avoid bosonization altogether as with the SCQRPA approach. 

Also in the work by Passos and al. [^8) the SCRPA method was applied to the present model. However, only the 
non superfluid formulation, i.e. SCRPA, was studied. The results are quite analogous to ours. In addition in |||] 
a further approximation, half way between RPA and SCRPA, the so-called renormalized RPA (r-RPA) where only 
the single-particle occupation numbers are allowed to be affected by ground-state correlation, has been considered. 
The astonishing finding there was that the exact occupation numbers are almost perfectly reproduced over the whole 
range of the coupling constant with r-RPA but not with SCRPA which under shoots the correlations. This was 
interpreted in |2^| as a positive feature of r-RPA over SCRPA. We can not follow this conclusion from our experience 
with SCRPA in this and other works ||). As we outlined above, any relax of the severe constraints of the Pauli 
principle respected in SCRPA will inevitably lead to more correlations as there should be. It can happen by accident 
that if one relaxes the Pauli principle just by the right amount that one falls more or less on the exact values. This 
is what happened for the QRPA ground-state energy discussed above and apparently it is also what happens for the 
occupation numbers in r-RPA. However, we think this result can not be generalized and for other models or physical 
situations the scenario may be completely different. The only really trust worthy theory is the full SCRPA approach, 
since it can be derived from a variational principle. If the results are not good, one must improve on SCRPA (i.e. 
include e.g. higher configurations) and not approximate it. 

In the work by Hagino and Bertsch |24| the QRPA approach is advocated. This in the spirit to have a numerically 
viable alternative to projected BCS and the method by Nogami-Lipkin J2{|. It is certainly true that in realistic cases 
SCRPA is numerically very demanding, though probably not impossible to solve with modern computers. Then, of 
course, in a first step it is worthwhile to investigate standard QRPA. This is for instance true if one intends to do 
large scale calculations for a great number of nuclei p4| . However, one should always remember that standard QRPA 
may have quite important failures which certainly will be most prominent in situations where the system is close to 
a phase transition. 



VI. THE QUESTION OF A SECOND CONSTRAINT ON THE PARTICLE NUMBER VARIANCE 

As we have seen above, with respect to BCS the SCQRPA reduces the spread in particle number by an important 
factor. However, the variance AA is still appreciable and one can ask the question whether it is not possible to 
further improve the theory on this point. A natural idea which comes to mind is that instead of fixing only (N) = N, 
one could at the same time fix (NN) = A^ 2 with a second Lagrange multiplier. Since in SCQRPA the number of 
variational parameters is largely increased with respect to BCS one could imagine that there is indeed enough freedom 
for constraining the particle number fluctuation to zero. The Hamiltonian to be considered is therefore 

H' = H - mN - n 2 N 2 . (67) 

Let us immediately give our conclusion : in the two level pairing case we could not find a solution to this problem. 
The system of non linear equations with the two constraints fj,i and Hi is quite complex and in spite of considerable 
numerical effort we did not have success to get the solution converged. We were not able to decide whether the 
difficulty is purely numerical or whether there is a principal problem. In fact we were at first encouraged by results 
we obtained in the one level pairing case (the seniority model). The outcome of employing the second constraint was 
that the one level model was solved exactly. In spite of being a somewhat trivial model which certainly limits the 
conclusions, it may be interesting to show how this goes. The Hamiltonian to be considered is now 



H' = H - mN - n 2 N 2 . 



(68) 
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where fix and \i 2 are two Lagrange multipliers fixing (N) = N and (N 2 ) = N 2 and 

H = -gQA^A (69) 

with in analogy to (^) = X) m >o a m a -mi anc ^ where we put the origin of energy at the single-particle level. As 
in the two level case we transform to quasi-particles and with only one level the SCQRPA equation reduce to a (2 x 2) 
eigenvalue problem 

{-b -x) (?)-*(?)= x2 " y2 = 1 (70) 

where in analogy with fl37|) 

r i 2 (#«> i ffi> i 

A = 2{g-A^ 2 )\{2v 2 {l-v 2 )-l)XY - v 2 {I - v 2 ){\ + 2Y 2 -Q a n' ) I (71) 



1 _ W 
1 n 



1 - 2 



B = 2(.g-4 M2 )<j {6v 2 (v 2 -1) + 1)XY -v 2 (l-v 2 )(l + 2Y 2 -il " 7" " 2 ) }. (72) 



1 n 



The Bogoliubov transformation to quasi-particles is obtained as for the case of two levels 



h A 

A -h I \v 



11 ' : e ( U ) (73) 



with in close analogy to the expression (|39|) of the two levels 

h = -(g + 4fi 2 n)v 2 - m, (74a) 

A = (ffO + 4^a)w - (g - A^)uv^2{XY + Y 2 ) + x ^ j, (74b) 

e = V^ 2 + A 2 . (74c) 



In addition to the SCQRPA equations (70) we have two further equations which, in principle, allow us to find the 



Lagrange multipliers [i\ and \i 2 (see, however, below) 

N = (TV) = (u 2 — v 2 )(N q ) + 2Qv 2 , (75a) 



N 2 = (N 2 ) = (u 2 -v 2 ) 2 (N 2 )+8nu 2 v 2 (l-^-){XY + Y 2 )+4nv 2 (u 2 + nv 2 ) 

+ 4v 2 (n{u 2 -v 2 ) -u 2 ){N q ) " (75b) 

We again see that eqs (75) reduce to the standard expressions, once, as in the HFB approximation, we pose Y = 



(N q ) = (N 2 ) — 0. In the case of the seniority model the number equation (75a) in the HFB approximation determines 



the amplitudes u, v and then no freedom is left to impose AN — y (N 2 ) — (N) — 0. However, in the more general 
approach of SCQRPA there is more freedom and, as we will see, we will be able to satisfy the relation AN = as 
well. For N q and N 2 we have the same relation as in ( |5l| ) and (^3|). Again the system of equations is therefore closed. 



Usually the number equation ( 75a ) and (75b) are to be used for the determination of the chemical potential fX\ and 
the second Lagrange multiplier [i 2 and the mean field equations for the amplitudes u, v and X, Y. In the present case 
it is, however, more convenient to invert the role of mean field and number equations, since eqs (|74|) do not depend 
on the Lagrange multipliers and therefore readily allow to determine v 2 and Y 2 as a function of the particle number 
N. Inversely the two mean field eqs (jfuj), ([73]) are linear in fXi and \i% and for instance it is seen that (f70|) directly 
yields 

V* = f (76) 
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independent of the particle number N. Considering, the well known exact expression for the ground-state energy of 
the model [§5) 

£ = -|(r> + l)iV+|iV 2 (77) 

we see from (X2 — d 9 ^§ that (f76| ) gives the exact value for the second Lagrange multiplier fj,2- For the chemical 
potential [i\ we obtain from ([73|) 



Mi = (.9 ~ 4 M2 )|(fi - l)v 2 + (1 - 2v 2 )(XY + Y 2 + - [ q) ^-Q ) I - | fi - 2/x 2 . (78) 

With relation ( f73| ) this gives /ii = — |(fi + 1) which again is the exact value. Furthermore, with ( |76| ) we have 
from (|l|), © that A = B = and therefore the RPA eigenvalue E = 0. This means that, as in standard QRPA, 
SCQRPA yields a Goldstone mode at zero energy.This feature is very rewarding, since it signifies that the particle 
number symmetry is exactly restored. 

It is well known that restoration of good particle number implies in this very simple model case that the model is solved 
exactly f25|| . We have already seen that one obtains the exact values for /ii and ^2- We now will show that one also 
obtains the exact value for the ground-state energy (and therefore for the whole band of ground-state energies). This 
goes as follows. For the expectation value of H of eq (|69| ) in the RPA ground-state, using the analogous relations ( |2g| ) 
and (B3) for this case and the quasi-particle representation for H, we can write 



E = (H) 

9 -(n + l){(l-2v 2 )(N q )+2nv 2 } + 9 



(l-2v 2 ) 2 (N) 



2 v- ■ -nv- - /,-'«/ ■ — > ' 4 

+ [Av 2 (n(l - 2v 2 ) -l + v 2 )- 8(1 - v 2 )v 2 {XY + Y 2 )](N q ) 
+ 80(1 - v 2 )v 2 (XY + Y 2 ) + 4Qv 2 (l - v 2 + nv 2 ) 



(79) 



In this expression we have used the Casimir relation for this case 4Y 2 (fl — (N q )) = 2(17 + l)(N q ) — (N 2 ) which follows 
from (El). Using the expression for TV and N 2 of (75a) and (75b) once more, we see that the exact expression d77 



recovered. It should be mentioned that because of the simplicity of the model also the Lipkin-Nogami approachEt 
solves the model exactly. 



VII. CONCLUSION 



In this work we extended for the first time the Self-Consistent RPA theory (SCQRPA) to the superfluid case for a 
model with more than one level. Indeed in [jl4| SCQRPA was already applied to the seniority Model but this only 
allowed to study rotation in gauge space whereas intrinsic excitations ("/3- vibrations") are absent in the + -sector of 
the seniority Model. We have considered the two level version of the pairing Hamiltonian with arbitrary degeneracies 
and fillings of the levels. We mostly considered the case J — 11/2 for the upper and lower levels. This configuration was 
chosen in order to have a better comparison with the work by Sambataro and Dinh Dang p3| ] which in many aspects 
is quite analogous to ours. Indeed SCQRPA can be considered as a Bogoliubov transformation among quasi-particle 
pair operators a^a^ and aa whereas in [^3| the quasi-particle pair operators were replaced by ideal bosons a^a^ ^ 
and then a Bogoliubov transformation among these boson operators was applied while the pairing Hamiltonian was 
also bosonized up to fourth order. For such collective pairs as they are formed in J = 11/2 shells a bosonization seems 
indeed valid and as expected our results are very close to the ones given in p3|] , even though they are consistently 
slightly better. This could be due to the fact that in SCQRPA one never bosonizes and rather all constraints from 
Pauli principle are fully kept. However, we do not want to attribute too much importance to these differences which 
only could become relevant for cases where a bosonization fails. On the other hand in our work considerably more 
issues were studied. In first place this concerns the physical interpretation and identification of the low lying state in 
SCQRPA. This state corresponds to the Goldstone mode in standard QRPA. However in SCQRPA this state comes at 
finite energy and reproduces very precisely the difference 2(/i + — /i~) of the chemical potentials of the exact solution. 
We also evaluated the fluctuation AN of the particle number and showed that with respect to the fluctuation in BCS 
theory there is a strong improvement. However, particle number symmetry is still not entirely restored. In spite of 
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this shortcoming for AN, for other quantities SCQRPA is always superior to BCS and QRPA approaches as explained 
in the main text. In fact the situation with respect to the particle number symmetry is somewhat particular and 
not encountered in other cases of spontaneously broken symmetries. For example in the case of rotation the angular 
momentum operator L z has no contributions which are hermitian in the deformed basis and then the Goldstone mode 
also comes in the case of SCRPA pjj . In order to improve on the restoration of particle number symmetry we also 
investigated the possibility of fixing (^V 2 ) = N 2 with a second Lagrange multiplier. Whereas in the one level pairing 
model we could show analytically that this solves the model exactly, in the two level model we could not find a 
numerical solution of the system of equations. It remained unclear whether this is due to some fundamental problem 
or just a numerical difficulty. 

We also discuss carefully in this work the transition from the non superfluid regime to the superfluid one. We for 
instance pointed out that the transition from SCRPA to SCQRPA is not continuous and in fact in both regimes quite 
different physical excitations are described. This also can be seen looking at the ground-state energies as a function 
of the coupling constant. In the transition region there is no continuous transition between the SCRPA and SCQRPA 
values but it is definitively seen that the SCRPA values for the ground-state energies deviate quite strongly from the 
exact values after the phase transition whereas SCQRPA stays close to them. 

In conclusion we can say that we have applied with very promising success for the first time SCQRPA to a more 
level pairing situation where, at least for the + -sector, all the complexity of a more realistic situation is present. It 
could be interesting to extend this work to the description of ultrasmall superconducting metallic grains for which 
the many-level picket fence model seems appropriate ||. 
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APPENDIX A: SOME USEFUL MATHEMATICAL RELATIONS 

At first, we will explain how we calculated the expectation values of type (Pj Pj'); we recall that j = ±1. From (|3^), 
i.e. the inversion of the excitation operators Q v , we can write, 

Pj = ^ Xj^Qv + Yj,uQ\,i 

v=l,2 

P i = E Xj,vQl+Yi,»Qu. (Al) 

!/=l,2 

Let us calculate for example (Pj Pj'), 



(P]Py) = V 1 " ng^f E Xj>AP}Q»)+Yi>AP}Qi)) (A2) 



' !/=l,2 

the first term in the rhs, is zero since Q v \) = 0. Therefore, we obtain, 



i>=1,2 



' ' v=l,2 

using the definition of the excitation operators (H1)j we find, 



E ^A[P],QU) (A3) 
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u=l,2 



= E^w-^X 1 -^)- ( A4 ) 

i/=l,2 ' 

We now will explain how we express the occupation number for each level j as function of the RPA amplitudes. We 
start with expectation value of (|5l]) in the RPA state, we can write 

(N q ,)^2(P]P J ) + I ^- i (pfpf). (A5) 

Using ph and <M) we find, 



where, 



(pfpf) = K jA + K 3 . 2 {P}P 3 ) + K^PjPj) + K M {PJP}) (A6) 



K i* = ^{-'(.V,.^.., + X ja Y ja ) 2 + :i(V;;, + V^fj. 

«i,4 = ^(A\:V,: • A,,V,. 2 )(Vr, • (A7) 



and, 



(Pi Pi) = (Y^+Y^l-^fl), 
(P}P}) = (X^+X^Y^il-^), 

(PjP}) = (A| 1+ A| 2 )(l-^1). 



Explicitly, we obtain 



K 4 (X jA Y jA +X j>2 Y j>2 ) \(1 - ^A). (A-) i 



Therefore, we can express (N q j) as function of the RPA amplitudes 



2(Y? 1+ Y? 2 )- 




+ ^ 2 ) 


+ ^ 3 (A| 1 




+ A 4 (x,, 1 y J , 1 






-^ 9 ) + ^{^i + 




-tf 3 (*j ! 1 i + 




- Ki(X jt xY u -f 





In the same way, we express (JV? and (iV 9i iA g _i) as follows, 



(iV^) = 2(0 + l)(A gj ) - 40(y^ + _ i^) (AH) 



(A8) 



(N q ,j) * ^ p (A10) 
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and. 



(JV,,iJV g ,_i) ~ 4(P 1 t P 1 Pl 1 P_ 1 



AM 



(A12) 



where M is a constant depending of the RPA amplitudes, it is given by 



M = 



(Yi 2 i + n 2 2)(^-i,i + ^-1,2) + (1 + ^) jV-i,i*i,i + Y_ l!2 X 1>2 ){X_ 1A Y 1A + X_i i2 F 1>2 ) + (^-1,1*1,1 
^-i,2^i, 2 )(X_ lil X 1 , 1 + X_ 1 , 2 X 1 , 2 ) j - ^(F-liXx,! +y_ 1 , 2 X 1)2 )(X_ 1 , 1 X 1 , 1 + X_ 1 , 2 X 1)2 )(P 1 t P 1 t ) 

(i r -i,iyi,i + y- 1)2 ii 1 2)(A'_i,iyi,i + x_i )2 yi,2)<PiPi) + + y_ 1 , 2 x 1 , 2 )(x_ 1 , 1 y 1 ,i 

X_i, 2 yi, 2 ) + (y-i.iFi.x +Y- 1 , 2 Y 1 , 2 )(X- 1 , 1 X 1 , 1 +X_i, 2 X 1 , 2 )|(2(P 1 t Pi) + (PPi 1 )) 



1 Ji -i,z i i.z; 1 — 1,1 1 J-i-L 1 1 -i,z i i,z/\-' l -i.i Jl i,i 1 - ,i -i^ ji i,z/ i j 1/ 1 \* 1 / j 

+ (Y-i.! + ^nl^-i^u + x_ 1 , 2 y 1 , 2 )(p 1 t Pl 1 ) + (y_ 2 M + y_\ 2 )(y_ M y M + y_ 1 , 2 y 1 , 2 )(p 1 t p_ 1 

+ \(Yli +y^ 2 )(y 2 1 ,i + y 2 i, 2 )(W, 1 ) + (aw)) 1 



(A13) 



APPENDIX B: METHOD FOR CALCULATION OF 

In this Appendix, we will present our method inspired from [^lj (see also p2|) for the derivation of the quantities 
of type N™j in the case of the two-levels pairing model. At first, we recall that in this model, the operators N q j, 

Pj ■ and P q j close the SU(2) algebra for each level j. Consequently the two level model fulfills an SU(2) x SU(2) 
algebra. Thanks to this special groupe structure, it is easy to find an complete ortho-normalized basis for the Hilbert 
subspace corresponding to each level j, 



1 Cl n i (SI - rij)\ p ^n s 



>, J = ±l 



where rij = 0, 1, . . . , SI. Using this basis, we can express the identity operator relatively to each level j as 



(Bl) 



<> 



n 



i;=£l»i>tel = l->H+£ 



7lj— n j— 1 

therefore, we can express the projector |— )(— | as follows 



!i LLpl 3 _\/_p™3 7 = ±1 

SV.nA 3 



n 



HH = ii-E n!n . 



One sees that (B3) produces an expansion of the form 



-)(-!= E ^Pprp., 



7Ua —0 



if we substitute (B4) in both lhs and rhs of (B3), we obtain the coefficients /3 m . 



0o = 1, An 3 = - E 



I. J. ; 



1=0 



si\(m 3 -i)\ 



(B2) 



(B3) 



(B4) 



(B5) 



For example, the first terms (3 mj are explicitly given by, 



n 2 -6^ + 12 

6(0- 2)' 



However, to calculate the quantities N q j and N%j, one can expand these operators as 

Q 

^-V/'pt 11 ^ 7 =±1 
l j = l 

For all operators of the form N k , using the fact that 



we can calculate 



N q,j\ n j) = 2n j\ n j)i 



N, 



= E^ fe >^ 

n i= 

= ^(2n J -) fc K-)(r 



rij — 

E 

Tl.- — 



{2n 3 fP] J \-)(-\P 



By inserting of (B4) in the rhs (B9c) and substituting (B7) into lhs of (B9c), we obtain 



N k - = V 

9,3 Z_^i 



Kj=0 



n^(fi-nj)! 
0!m! 



rrij—O 



{2{l j -m j )fp mj pf 1 P]- 



l j= l n,=l """ ■' 

- EE 



0!n,! 

m_j — 

n l s- m i{Q,-lj+mj)\ 



ij = l ITlj=0 



Therefore, by identification, from ( |B10| ) it is easy to get the coefficients a[ k ^ , 



z, — 1 



x Qh-">-i (fi - lj +mj)\ 



To calculate iV gj - we put k = 1 in (B7), and find 



with, 



(i) _ f^" TO * (0 - lj +mj)\ 



E 

mi —0 



QlQj — TTij ) ! 



(2(ij - mjJJA, 
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The first coefficients af^ are explicitly given by, 

a? = 2, 
a? 



2 



o- 1' 

4 



d (0-l)(0-2)' 
(1) 2(50-6) 

4 (0-l) 2 (0-2)(0-3)' 1 ; 



Therefore, for N q j we can write 



+ 2 + 2 9 4 +3 , 2(50-6) +4 - 

J . = opfp i pT p2 j pt p3 i V )_ pt p4 

,,J J o — i J 1 (n-i)(n-2) J J (o-i) 2 (o-2)(o-3) J J 



In the same way, it is very easy with this method to find an expansion for JVf •, it is sufficient to put /e = 2 in (B7) 



f 2 1 ) 

and calculate a) . We find 



+ 4(0 + 1) +2 „ 8(0 + 1) + 3 „ 4(0 + l)(50-6) +4 a , 

N ■— 4P T P • -I ' — -P! P A - — ' — P T P 3 -I ' — — P T P 4- (B16) 

« = J+ 0-1 ^ J + (0-l)(0-2) * i + (0-l)2(0-2)(0-3) * * 1 j 
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FIG. 1: Ground-state energy E ga and excitation energy of the first + state E exc as a function of the variable V = gQ,/2e 
described in the text and for particle number N = 12 (energies are divided by 2e). The spin of the levels is J — 11/2. The 
results refer to exact calculations (solid line and double-dot dashed line), BCS (dotted line), RPA and QRPA (dot-dashed line), 
SCQRPA (dashed line) and SCRPA (double-dash dotted line). (Note that SCRPA and SCQRPA solutions co-exist over a wide 
range of V- values). 
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FIG. 2: Zoom on the ground-state energy E gs as a function of the variable V = gQ/2e described in the text and for particle 
number N = 12 (energies are divided by 2e). The spin of the levels is J = 11/2. The results refer to exact calculations (solid 
line), BCS (dotted line), RPA and QRPA (dot-dashed line), SCQRPA (dashed line) and SCRPA (double-dash dotted line). 
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FIG. 4: As in Fig.| but for N = 4. 
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FIG. 6: Excitation energy of the soft (spurious) mode (energies are divided by 2e) as a function of the variable V = gQ./2e 
described in the text and for particle number N = 10. The spin of the levels is J = 11/2. The results refer to exact calculations 
(solid line, dotted line and double-dot dashed line), QRPA (dot-dashed line), and SCQRPA (dashed line). 





1 1 1 1 

N+i N-2 N 

— p 1N + P - 2 F 
p N+2 p N 


1 1 1 

N = 4 






p N ~ 2 p N 








SCQRPA 








_ QRPA 








1 , 1 


1 , 1 





5 l i l i l i l i l_ 

0.2 0.4 0.6 0.8 

V 

FIG. 7: As in Fig.| but for N = 4. 
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FIG. 8: Excitation energy of the soft (spurious) mode (energies are divided by 2e) as a function of the variable V = gO,/2e 
described in the text and for particle number iV = 12. The spin of the levels is J = 11/2. The results refer to exact calculations 
(solid line, dotted line and double-dot dashed line), RPA and QRPA (dot-dashed line), SCQRPA (dashed line) and SCRPA 
(double-dash dotted line). 




FIG. 9: Comparison between SCQRPA, BCS and exact results for the chemical potentials /i = + fi— ) and fj, = 

±±(£^ ±2 - Eg), for particle number N = 4. The spin of the levels is J = 11/2. The results refer to exact calculations (solid 
line), SCQRPA (dashed line) and BCS (dotted line). 




FIG. 10: As in Fig.| but for N = 8. 




FIG. 11: As in Fig.| but for N = 12. 
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FIG. 12: Excitation energies 2/i + (upper lines) and 2fi~ (energies are divided by 2e) as a function of the variable V = gil/2e 
described in the text and for N = 12. The spin of the levels is J = 11/2. The full lines correspond to the exact results, the 
broken lines to SCRPA with occupation numbers calculated with the wave function ([io|), and dotted lines to SCRPA with 
occupation numbers from the expansion (Ell) . 




FIG. 13: Particle number in the upper level iVi as function of the variable V = gfl/2e described in the text and for = 12. The 
spin of the levels is J = 11/2. The results refer to exact calculations (solid line), SCQRPA (dashed line), SCRPA (double-dash 
dotted line) and QRPA (dot-dashed line). 
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FIG. 14: Variance as a function of the variable V = gQ./2e described in the text and for particle number N = 10. The spin of 
the levels is J = 11/2. The results refer to SCQRPA calculations (dashed line) and BCS (dotted line). 
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FIG. 15: As in Fig|| but N = 12. 
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FIG. 17: Correlation lunction C as a function of the variable V = gQ,/2e described in the text and for particle number N — 12. 
The spin of the levels is J = 11/2. The results refer to SCQRPA (solid line), renormalized BCS (dashed line) and standard BCS 
(dotted line). 



